#
# Created 24 Apr 2019
#
# Replication file for:
## Clark, Zucker, and Urpelainen.
## "Political Institutions and Pollution: Evidence from Coal-Fired Power Generation."
## Review of Policy Research.
#

########## #
# SETUP ####
########## #

rm(list=ls())

pacman::p_load(dplyr,
               foreign,
               ggplot2,
               ggmap,
               interplot,
               lattice,
               lmtest,
               mapdata,
               maps,
               multiwayvcov,
               nlme,
               plm,
               raster,
               rworldmap,
               stargazer,
               xtable)

# > Data loading ####

coal <- read.csv("rep_data.csv")

############### #
# MAIN PAPER ####
############### #

# > Figure 1: Map ####

wmdat <- coal[coal$year == 2016,]
wmdat <- wmdat[,c("year","country","population","cumul_capacity_post80")]
wmdat$country <- as.character(wmdat$country)

map.world <- map_data("world")

wmdat$country[wmdat$country == "United States"] <- "USA"
wmdat$country[wmdat$country == "United Kingdom"] <- "UK"

map.world_joined <- left_join(map.world, wmdat, by = c("region" = "country"))
map.world_joined$cumul_capacity_post80[is.na(map.world_joined$cumul_capacity_post80)] <- 0

map.world_joined$cumul_capacity_post80_LOG <- log1p(map.world_joined$cumul_capacity_post80)

map.world_joined$cumul_capacity_post80_PC <- map.world_joined$cumul_capacity_post80/map.world_joined$population

map.world_joined$cumul_capacity_post80_PC[is.na(map.world_joined$cumul_capacity_post80_PC)] <- 0

map_output <- ggplot() +
  geom_polygon(data = map.world_joined, aes(x = long, y = lat, group = group, fill = cumul_capacity_post80_LOG), color = "black") +
  scale_fill_gradient(low = "white", high = "red1", name = "Megawatts (log)"
  ) +
  theme(text = element_text(color = "black")
        ,panel.background = element_rect(fill = "white")
        ,plot.background = element_rect(fill = "white")
        ,legend.background = element_rect(fill = "white")
        ,panel.grid = element_blank()
        ,plot.title = element_text(size = 20)
        ,plot.subtitle = element_text(size = 10)
        ,axis.text = element_blank()
        ,axis.title = element_blank()
        ,axis.ticks = element_blank()
        ,legend.position = "right"
  )
map_output

# > Figure 2: Plot ####

cumul_cap <- coal[,c("year","cumul_capacity_post80","cgv_democracy","population")]
cumul_cap <- aggregate(cbind(cumul_capacity_post80, population) ~ cgv_democracy + year, data = cumul_cap, FUN = sum) #NOTE: only summing cumulative capacities for country-years where population data is available
cumul_cap$cgv_democracy <- as.character(cumul_cap$cgv_democracy)
cumul_cap$cumul_pc <- (cumul_cap$cumul_capacity_post80-1)/cumul_cap$population

cumul_cap2 <- coal[coal$country != "China",]
cumul_cap2 <- cumul_cap2[,c("year","cumul_capacity_post80","bmr_democracy","population")]
cumul_cap2$cumul_capacity_post80 <- cumul_cap2$cumul_capacity_post80 - 1
cumul_cap2 <- aggregate(cbind(cumul_capacity_post80, population) ~ bmr_democracy + year, data = cumul_cap2, FUN = sum) #NOTE: only summing cumulative capacities for country-years where population data is available
cumul_cap2$bmr_democracy <- as.character(cumul_cap2$bmr_democracy)
cumul_cap2$cumul_pc <- ((cumul_cap2$cumul_capacity_post80)/cumul_cap2$population) * 1000000

cumul_cap2_china <- coal[coal$country == "China",]
cumul_cap2_china <- cumul_cap2_china[,c("year","cumul_capacity_post80","bmr_democracy","population")]
cumul_cap2_china$cumul_capacity_post80 <- cumul_cap2_china$cumul_capacity_post80 - 1
cumul_cap2_china <- aggregate(cbind(cumul_capacity_post80, population) ~ bmr_democracy + year, data = cumul_cap2_china, FUN = sum) #NOTE: only summing cumulative capacities for country-years where population data is available
cumul_cap2_china$bmr_democracy <- as.character(cumul_cap2_china$bmr_democracy)
cumul_cap2_china$cumul_pc <- ((cumul_cap2_china$cumul_capacity_post80)/cumul_cap2_china$population) * 1000000
cumul_cap2_china$bmr_democracy <- 2

cumul_cap2 <- rbind(cumul_cap2, cumul_cap2_china)
cumul_cap2$cumul_ln_pc <- log(cumul_cap2$cumul_capacity_post80/cumul_cap2$population * 1000000)

cumul_plot2 <- ggplot(cumul_cap2, aes(x = year, y = cumul_pc, group = bmr_democracy, color = bmr_democracy)) +
  geom_line()
cumul_plot2 <- cumul_plot2 + scale_color_brewer(palette="Set1", name = "Legend", labels = c("Non-Democracy","Democracy","China"))
cumul_plot2 <- cumul_plot2 + theme_minimal() +
  xlab("Year") + ylab("Megawatts per millions of people")
cumul_plot2

# > Figure 3: Plot ####

cumul_plot3 <- ggplot(cumul_cap2, aes(x = year, y = cumul_ln_pc, group = bmr_democracy, color = bmr_democracy)) +
  geom_line()
cumul_plot3 <- cumul_plot3 + scale_color_brewer(palette="Set1", name = "Legend", labels = c("Non-Democracy","Democracy","China"))
cumul_plot3 <- cumul_plot3 + theme_minimal() +
  xlab("Year") + ylab("Megawatts per millions of people (ln)")
cumul_plot3

# > Figure 4: Interaction Plot ####

reg2_int <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5 +
                     log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                     urban_pop_l5 +
                     kyoto_l5 + EU_l5,
                   data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int <- coeftest(reg2_int, vcov=function(x) vcovHC(x, cluster = "group")) )

int1.vcmat <- vcovHC(reg2_int, cluster = "group")
int1.covariance <- int1.vcmat["polity2_l5","polity2_l5:oeec_oecd_l5"]
int1.var1 <- int1.vcmat["polity2_l5","polity2_l5"]
int1.var2 <- int1.vcmat["oeec_oecd_l5","oeec_oecd_l5"]
int1.var3 <- int1.vcmat["polity2_l5:oeec_oecd_l5","polity2_l5:oeec_oecd_l5"]

int1.df <- data.frame(fake = rep(NA,2), coef1 = rep(NA,2), lb = rep(NA,2), ub = rep(NA,2))
int1.df$fake <- 0:1 # values of OECD membership
int1.df$coef1 <- reg2_int$coefficients[["polity2_l5"]] + 
  reg2_int$coefficients[["polity2_l5:oeec_oecd_l5"]]*int1.df$fake # coefficients for each OECD value
int1.df$lb <- int1.df$coef1 - 1.96*sqrt(int1.var1 + (int1.df$fake)^2*int1.var3 + 2*int1.df$fake*int1.covariance) # lower bound of 95% CI
int1.df$ub <- int1.df$coef1 + 1.96*sqrt(int1.var1 + (int1.df$fake)^2*int1.var3 + 2*int1.df$fake*int1.covariance) # upper bound of 95% CI

int1 <- interplot(int1.df, ercolor = "gray", esize = 1.5) +
  geom_point(size = 2, color = "red") + ylab("Estimated Coefficient of Polity2 (t-5)") +
  geom_line(linetype = "solid", color = "red") + xlab("OECD Membership (t-5)") + theme_bw() + 
  theme(legend.position="none") + geom_hline(yintercept = 0, linetype = "dashed")
int1

# > Figure 5: Interaction Plot ####

coal2 <- coal
coal2$log_gdp_pc_l5 <- log(coal2$gdp_pc_l5)

reg2_int2_mfx <- plm(log(cumul_capacity_post80) ~ polity2_l5*log_gdp_pc_l5 +
                          log(population_l5) + trade_pctgdp_l5 +
                          urban_pop_l5 +
                          kyoto_l5 + EU_l5 + oeec_oecd_l5,
                        data = coal2, model = "within", effect = "twoways", index = c("country", "year"))

int2.vcmat <- vcovHC(reg2_int2_mfx, cluster = "group")
int2.covariance <- int2.vcmat["polity2_l5","polity2_l5:log_gdp_pc_l5"]
int2.var1 <- int2.vcmat["polity2_l5","polity2_l5"]
int2.var2 <- int2.vcmat["log_gdp_pc_l5","log_gdp_pc_l5"]
int2.var3 <- int2.vcmat["polity2_l5:log_gdp_pc_l5","polity2_l5:log_gdp_pc_l5"]

int2.df <- data.frame(fake = rep(NA,1000), coef1 = rep(NA,1000), lb = rep(NA,1000), ub = rep(NA,1000))
int2.df$fake <- seq(from = 5.5, to = 11, by = (11-5.5)/999) # values of log(GDP per capita)
int2.df$coef1 <- reg2_int2_mfx$coefficients[["polity2_l5"]] + 
  reg2_int2_mfx$coefficients[["polity2_l5:log_gdp_pc_l5"]]*int2.df$fake # coefficients for each log(GDP per capita) value
int2.df$lb <- int2.df$coef1 - 1.96*sqrt(int2.var1 + int2.df$fake^2*int2.var3 + 2*int2.df$fake*int2.covariance) # lower bound of 95% CI
int2.df$ub <- int2.df$coef1 + 1.96*sqrt(int2.var1 + int2.df$fake^2*int2.var3 + 2*int2.df$fake*int2.covariance) # upper bound of 95% CI

gdppc.hist <- coal2$log_gdp_pc_l5[complete.cases(coal2$cumul_capacity_post80,
                                                 coal2$polity2_l5, coal2$oeec_oecd_l5,
                                                 coal2$gdp_pc_l5, coal2$population_l5,
                                                 coal2$trade_pctgdp_l5, coal2$urban_pop_l5,
                                                 coal2$kyoto_l5, coal2$EU_l5)]

int2 <- interplot(int2.df, var2_dt = gdppc.hist, hist = TRUE)  +
  aes(colour = "blue") + xlab("GDP/capita (ln, t-5)") +
  ylab("Estimated Coefficient of Polity2 (t-5)") + theme_bw() + 
  theme(legend.position="none") + geom_hline(yintercept = 0, linetype = "dashed")
int2

# > Table 1: Summary Statistics ####

coal_analysis <- coal[coal$year >= 1980 & coal$year < 2017,]

coal_sumstats <- coal_analysis[,c("year","cumul_capacity_post80","polity2_l5",
                          "oeec_oecd_l5","gdp_pc_l5","population_l5","urban_pop_l5","trade_pctgdp_l5",
                          "kyoto_l5","EU_l5")]
coal_sumstats_vn <- c("Year","Cumulative commissioned capacity",
              "Polity2 ($t-5$)", "OECD membership ($t-5$)", "GDP/capita ($t-5$)",
              "Population ($t-5$)", "Urbanization, \\%  ($t-5$)","Trade, \\% GDP ($t-5$)",
              "Kyoto ratification ($t-5$)", "EU membership ($t-5$)")

stargazer(title = "Summary Statistics",
          coal_sumstats,
          header = F,
          covariate.labels = coal_sumstats_vn,
          omit.summary.stat = c("p25", "p75"),
          float = F)

# > Table 2: Main Regression Table ####

reg1 <- plm(log(cumul_capacity_post80) ~ polity2_l5,
               data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1 <- coeftest(reg1, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2 <- plm(log(cumul_capacity_post80) ~ polity2_l5 +
                 log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                 urban_pop_l5 +
                 kyoto_l5 + EU_l5 + oeec_oecd_l5,
               data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2 <- coeftest(reg2, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5,
                   data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int <- coeftest(reg1_int, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5 +
                     log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                     urban_pop_l5 +
                     kyoto_l5 + EU_l5,
                   data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int <- coeftest(reg2_int, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2 <- plm(log(cumul_capacity_post80) ~ polity2_l5*log(gdp_pc_l5),
                    data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2 <- coeftest(reg1_int2, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2 <- plm(log(cumul_capacity_post80) ~ polity2_l5*log(gdp_pc_l5) +
                      log(population_l5) + trade_pctgdp_l5 +
                      urban_pop_l5 +
                      kyoto_l5 + EU_l5 + oeec_oecd_l5,
                    data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2 <- coeftest(reg2_int2, vcov=function(x) vcovHC(x, cluster = "group")) )

list_main <- list(reg1,
                  reg2,
                  reg1_int,
                  reg2_int,
                  reg1_int2,
                  reg2_int2)
se_main <- list(coef_reg1[,2],
                coef_reg2[,2],
                coef_reg1_int[,2],
                coef_reg2_int[,2],
                coef_reg1_int2[,2],
                coef_reg2_int2[,2])
p_main <- list(coef_reg1[,4],
               coef_reg2[,4],
               coef_reg1_int[,4],
               coef_reg2_int[,4],
               coef_reg1_int2[,4],
               coef_reg2_int2[,4])
covlabs_main <- c("Polity2 ($t-5$)",
                  "GDP/capita, \\$ (ln, $t-5$)",
                  "Population (ln, $t-5$)",
                  "Trade, \\% GDP ($t-5$)",
                  "Urban population, \\% ($t-5$)",
                  "Kyoto ratification ($t-5$)",
                  "EU member ($t-5$)",
                  "OECD member ($t-5$)",
                  "Polity2 $\\times$ OECD member ($t-5$)",
                  "Polity2 $\\times$ GDP/capita (ln) ($t-5$)")
n_main <- list(c("Countries","71","69","71","69","69","69"))

omit.stats <- c("rsq","ser","f","adj.rsq")

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_main,
          header = F,
          list_main,
          covariate.labels = covlabs_main,
          omit.stat = omit.stats,
          add.lines = n_main,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

############### #
# APPENDICES ####
############### #

# > Cumulative coal capacity tables ####

data2016 <- coal[coal$year == 2016,]
data2016 <- as_data_frame(data2016)
data2016 <- data2016[,c("country","cumul_capacity_post80")]
data2016$cumul_capacity_post80 <- data2016$cumul_capacity_post80 - 1
data2016$cumul_capacity_post80 <- round(data2016$cumul_capacity_post80)
colnames(data2016) <- c("Country","Cumul. Capacity")

print(xtable(data2016[c(1:35),], digits = c(0,0,0)),
      floating=FALSE,
      font.size="footnotesize",
      latex.environments="center",
      include.rownames=F)

print(xtable(data2016[c(36:71),], digits = c(0,0,0)),
      floating=FALSE,
      font.size="footnotesize",
      latex.environments="center",
      include.rownames=F)

# > Three-year lag ####

reg1_int_t3 <- plm(log(cumul_capacity_post80) ~ polity2_l3*oeec_oecd_l3,
                      data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_t3 <- coeftest(reg1_int_t3, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_t3 <- plm(log(cumul_capacity_post80) ~ polity2_l3*oeec_oecd_l3 +
                        log(gdp_pc_l3) + log(population_l3) + trade_pctgdp_l3 +
                        urban_pop_l3 +
                        kyoto_l3 + EU_l3,
                      data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_t3 <- coeftest(reg2_int_t3, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_t3 <- plm(log(cumul_capacity_post80) ~ polity2_l3*log(gdp_pc_l3),
                       data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_t3 <- coeftest(reg1_int2_t3, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_t3 <- plm(log(cumul_capacity_post80) ~ polity2_l3*log(gdp_pc_l3) +
                         log(population_l3) + trade_pctgdp_l3 +
                         urban_pop_l3 +
                         kyoto_l3 + EU_l3 + oeec_oecd_l3,
                       data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_t3 <- coeftest(reg2_int2_t3, vcov=function(x) vcovHC(x, cluster = "group")) )

list_t3 <- list(reg1_int_t3,
                reg2_int_t3,
                reg1_int2_t3,
                reg2_int2_t3)
se_t3 <- list(coef_reg1_int_t3[,2],
              coef_reg2_int_t3[,2],
              coef_reg1_int2_t3[,2],
              coef_reg2_int2_t3[,2])
p_t3 <- list(coef_reg1_int_t3[,4],
             coef_reg2_int_t3[,4],
             coef_reg1_int2_t3[,4],
             coef_reg2_int2_t3[,4])
covlabs_t3 <- c("Polity2 ($t-3$)",
                "OECD member ($t-3$)",
                "GDP/capita, \\$ (ln, $t-3$)",
                "Population (ln, $t-3$)",
                "Trade, \\% GDP ($t-3$)",
                "Urban population, \\% ($t-3$)",
                "Kyoto ratification ($t-3$)",
                "EU member ($t-3$)",
                "Polity2 $\\times$ OECD member ($t-3$)",
                "Polity2 $\\times$ GDP/capita (ln) ($t-3$)")
n_t3 <- list(c("Countries","70","69","69","69"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_t3,
          header = F,
          list_t3,
          covariate.labels = covlabs_t3,
          omit.stat = omit.stats,
          add.lines = n_t3,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > Seven-year lag ####

reg1_int_t7 <- plm(log(cumul_capacity_post80) ~ polity2_l7*oeec_oecd_l7,
                      data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_t7 <- coeftest(reg1_int_t7, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_t7 <- plm(log(cumul_capacity_post80) ~ polity2_l7*oeec_oecd_l7 +
                        log(gdp_pc_l7) + log(population_l7) + trade_pctgdp_l7 +
                        urban_pop_l7 +
                        kyoto_l7 + EU_l7,
                      data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_t7 <- coeftest(reg2_int_t7, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_t7 <- plm(log(cumul_capacity_post80) ~ polity2_l7*log(gdp_pc_l7),
                       data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_t7 <- coeftest(reg1_int2_t7, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_t7 <- plm(log(cumul_capacity_post80) ~ polity2_l7*log(gdp_pc_l7) +
                         log(population_l7) + trade_pctgdp_l7 +
                         urban_pop_l7 +
                         kyoto_l7 + EU_l7 + oeec_oecd_l7,
                       data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_t7 <- coeftest(reg2_int2_t7, vcov=function(x) vcovHC(x, cluster = "group")) )

list_t7 <- list(reg1_int_t7,
                reg2_int_t7,
                reg1_int2_t7,
                reg2_int2_t7)
se_t7 <- list(coef_reg1_int_t7[,2],
              coef_reg2_int_t7[,2],
              coef_reg1_int2_t7[,2],
              coef_reg2_int2_t7[,2])
p_t7 <- list(coef_reg1_int_t7[,4],
             coef_reg2_int_t7[,4],
             coef_reg1_int2_t7[,4],
             coef_reg2_int2_t7[,4])
covlabs_t7 <- c("Polity2 ($t-7$)",
                "OECD member ($t-7$)",
                "GDP/capita, \\$ (ln, $t-7$)",
                "Population (ln, $t-7$)",
                "Trade, \\% GDP ($t-7$)",
                "Urban population, \\% ($t-7$)",
                "Kyoto ratification ($t-7$)",
                "EU member ($t-7$)",
                "Polity2 $\\times$ OECD member ($t-7$)",
                "Polity2 $\\times$ GDP/capita (ln) ($t-7$)")
n_t7 <- list(c("Countries","71","69","69","69"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_t7,
          header = F,
          list_t7,
          covariate.labels = covlabs_t7,
          omit.stat = omit.stats,
          add.lines = n_t7,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > EXREC regression ####

reg1_int_exrec <- plm(log(cumul_capacity_post80) ~ exrec*oeec_oecd_l5,
                         data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_exrec <- coeftest(reg1_int_exrec, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_exrec <- plm(log(cumul_capacity_post80) ~ exrec*oeec_oecd_l5 +
                           log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                           urban_pop_l5 +
                           kyoto_l5 + EU_l5,
                         data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_exrec <- coeftest(reg2_int_exrec, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_exrec <- plm(log(cumul_capacity_post80) ~ exrec*log(gdp_pc_l5),
                          data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_exrec <- coeftest(reg1_int2_exrec, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_exrec <- plm(log(cumul_capacity_post80) ~ exrec*log(gdp_pc_l5) +
                            log(population_l5) + trade_pctgdp_l5 +
                            urban_pop_l5 +
                            kyoto_l5 + EU_l5 + oeec_oecd_l5,
                          data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_exrec <- coeftest(reg2_int2_exrec, vcov=function(x) vcovHC(x, cluster = "group")) )

list_exrec <- list(reg1_int_exrec,
                   reg2_int_exrec,
                   reg1_int2_exrec,
                   reg2_int2_exrec)
se_exrec <- list(coef_reg1_int_exrec[,2],
                 coef_reg2_int_exrec[,2],
                 coef_reg1_int2_exrec[,2],
                 coef_reg2_int2_exrec[,2])
p_exrec <- list(coef_reg1_int_exrec[,4],
                coef_reg2_int_exrec[,4],
                coef_reg1_int2_exrec[,4],
                coef_reg2_int2_exrec[,4])
covlabs_exrec <- c("EXREC ($t-5$)",
                   "OECD member ($t-5$)",
                   "GDP/capita, \\$ (ln, $t-5$)",
                   "Population (ln, $t-5$)",
                   "Trade, \\% GDP ($t-5$)",
                   "Urban population, \\% ($t-5$)",
                   "Kyoto ratification ($t-5$)",
                   "EU member ($t-5$)",
                   "EXREC $\\times$ OECD member ($t-5$)",
                   "EXREC $\\times$ GDP/capita (ln) ($t-5$)")
n_exrec <- list(c("Countries","70","68","68","68"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_exrec,
          header = F,
          list_exrec,
          covariate.labels = covlabs_exrec,
          omit.stat = omit.stats,
          add.lines = n_exrec,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > EXCONST regression ####

reg1_int_exconst <- plm(log(cumul_capacity_post80) ~ exconst*oeec_oecd_l5,
                           data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_exconst <- coeftest(reg1_int_exconst, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_exconst <- plm(log(cumul_capacity_post80) ~ exconst*oeec_oecd_l5 +
                             log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                             urban_pop_l5 +
                             kyoto_l5 + EU_l5,
                           data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_exconst <- coeftest(reg2_int_exconst, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_exconst <- plm(log(cumul_capacity_post80) ~ exconst*log(gdp_pc_l5),
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_exconst <- coeftest(reg1_int2_exconst, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_exconst <- plm(log(cumul_capacity_post80) ~ exconst*log(gdp_pc_l5) +
                              log(population_l5) + trade_pctgdp_l5 +
                              urban_pop_l5 +
                              kyoto_l5 + EU_l5 + oeec_oecd_l5,
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_exconst <- coeftest(reg2_int2_exconst, vcov=function(x) vcovHC(x, cluster = "group")) )

list_exconst <- list(reg1_int_exconst,
                     reg2_int_exconst,
                     reg1_int2_exconst,
                     reg2_int2_exconst)
se_exconst <- list(coef_reg1_int_exconst[,2],
                   coef_reg2_int_exconst[,2],
                   coef_reg1_int2_exconst[,2],
                   coef_reg2_int2_exconst[,2])
p_exconst <- list(coef_reg1_int_exconst[,4],
                  coef_reg2_int_exconst[,4],
                  coef_reg1_int2_exconst[,4],
                  coef_reg2_int2_exconst[,4])
covlabs_exconst <- c("EXCONST ($t-5$)",
                     "OECD member ($t-5$)",
                     "GDP/capita, \\$ (ln, $t-5$)",
                     "Population (ln, $t-5$)",
                     "Trade, \\% GDP ($t-5$)",
                     "Urban population, \\% ($t-5$)",
                     "Kyoto ratification ($t-5$)",
                     "EU member ($t-5$)",
                     "EXCONST $\\times$ OECD member ($t-5$)",
                     "EXCONST $\\times$ GDP/capita (ln) ($t-5$)")
n_exconst <- list(c("Countries","70","68","68","68"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_exconst,
          header = F,
          list_exconst,
          covariate.labels = covlabs_exconst,
          omit.stat = omit.stats,
          add.lines = n_exconst,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > POLCOMP regression ####

reg1_int_polcomp <- plm(log(cumul_capacity_post80) ~ polcomp*oeec_oecd_l5,
                           data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_polcomp <- coeftest(reg1_int_polcomp, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_polcomp <- plm(log(cumul_capacity_post80) ~ polcomp*oeec_oecd_l5 +
                             log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                             urban_pop_l5 +
                             kyoto_l5 + EU_l5,
                           data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_polcomp <- coeftest(reg2_int_polcomp, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_polcomp <- plm(log(cumul_capacity_post80) ~ polcomp*log(gdp_pc_l5),
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_polcomp <- coeftest(reg1_int2_polcomp, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_polcomp <- plm(log(cumul_capacity_post80) ~ polcomp*log(gdp_pc_l5) +
                              log(population_l5) + trade_pctgdp_l5 +
                              urban_pop_l5 +
                              kyoto_l5 + EU_l5 + oeec_oecd_l5,
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_polcomp <- coeftest(reg2_int2_polcomp, vcov=function(x) vcovHC(x, cluster = "group")) )

list_polcomp <- list(reg1_int_polcomp,
                     reg2_int_polcomp,
                     reg1_int2_polcomp,
                     reg2_int2_polcomp)
se_polcomp <- list(coef_reg1_int_polcomp[,2],
                   coef_reg2_int_polcomp[,2],
                   coef_reg1_int2_polcomp[,2],
                   coef_reg2_int2_polcomp[,2])
p_polcomp <- list(coef_reg1_int_polcomp[,4],
                  coef_reg2_int_polcomp[,4],
                  coef_reg1_int2_polcomp[,4],
                  coef_reg2_int2_polcomp[,4])
covlabs_polcomp <- c("POLCOMP ($t-5$)",
                     "OECD member ($t-5$)",
                     "GDP/capita, \\$ (ln, $t-5$)",
                     "Population (ln, $t-5$)",
                     "Trade, \\% GDP ($t-5$)",
                     "Urban population, \\% ($t-5$)",
                     "Kyoto ratification ($t-5$)",
                     "EU member ($t-5$)",
                     "POLCOMP $\\times$ OECD member ($t-5$)",
                     "POLCOMP $\\times$ GDP/capita (ln) ($t-5$)")
n_polcomp <- list(c("Countries","70","68","68","68"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_polcomp,
          header = F,
          list_polcomp,
          covariate.labels = covlabs_polcomp,
          omit.stat = omit.stats,
          add.lines = n_polcomp,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > Boix et al. binary regression ####

reg1_int_bmr <- plm(log(cumul_capacity_post80) ~ bmr_democracy_l5*oeec_oecd_l5,
                       data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_bmr <- coeftest(reg1_int_bmr, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_bmr <- plm(log(cumul_capacity_post80) ~ bmr_democracy_l5*oeec_oecd_l5 +
                         log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                         urban_pop_l5 +
                         kyoto_l5 + EU_l5,
                       data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_bmr <- coeftest(reg2_int_bmr, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_bmr <- plm(log(cumul_capacity_post80) ~ bmr_democracy_l5*log(gdp_pc_l5),
                        data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_bmr <- coeftest(reg1_int2_bmr, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_bmr <- plm(log(cumul_capacity_post80) ~ bmr_democracy_l5*log(gdp_pc_l5) +
                          log(population_l5) + trade_pctgdp_l5 +
                          urban_pop_l5 +
                          kyoto_l5 + EU_l5 + oeec_oecd_l5,
                        data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_bmr <- coeftest(reg2_int2_bmr, vcov=function(x) vcovHC(x, cluster = "group")) )

list_bmr <- list(reg1_int_bmr,
                 reg2_int_bmr,
                 reg1_int2_bmr,
                 reg2_int2_bmr)
se_bmr <- list(coef_reg1_int_bmr[,2],
               coef_reg2_int_bmr[,2],
               coef_reg1_int2_bmr[,2],
               coef_reg2_int2_bmr[,2])
p_bmr <- list(coef_reg1_int_bmr[,4],
              coef_reg2_int_bmr[,4],
              coef_reg1_int2_bmr[,4],
              coef_reg2_int2_bmr[,4])
covlabs_bmr <- c("BMR democracy ($t-5$)",
                 "OECD member ($t-5$)",
                 "GDP/capita, \\$ (ln, $t-5$)",
                 "Population (ln, $t-5$)",
                 "Trade, \\% GDP ($t-5$)",
                 "Urban population, \\% ($t-5$)",
                 "Kyoto ratification ($t-5$)",
                 "EU member ($t-5$)",
                 "Polity2 $\\times$ OECD member ($t-5$)",
                 "Polity2 $\\times$ GDP/capita (ln) ($t-5$)")
n_bmr <- list(c("Countries","71","70","70","70"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_bmr,
          header = F,
          list_bmr,
          covariate.labels = covlabs_bmr,
          omit.stat = omit.stats,
          add.lines = n_bmr,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > Unified Democracy Scores regression ####

reg1_int_udsmean <- plm(log(cumul_capacity_post80) ~ uds_mean_l5*oeec_oecd_l5,
                           data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_udsmean <- coeftest(reg1_int_udsmean, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_udsmean <- plm(log(cumul_capacity_post80) ~ uds_mean_l5*oeec_oecd_l5 +
                             log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                             urban_pop_l5 +
                             kyoto_l5 + EU_l5,
                           data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_udsmean <- coeftest(reg2_int_udsmean, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_udsmean <- plm(log(cumul_capacity_post80) ~ uds_mean_l5*log(gdp_pc_l5),
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_udsmean <- coeftest(reg1_int2_udsmean, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_udsmean <- plm(log(cumul_capacity_post80) ~ uds_mean_l5*log(gdp_pc_l5) +
                              log(population_l5) + trade_pctgdp_l5 +
                              urban_pop_l5 +
                              kyoto_l5 + EU_l5 + oeec_oecd_l5,
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_udsmean <- coeftest(reg2_int2_udsmean, vcov=function(x) vcovHC(x, cluster = "group")) )

list_uds <- list(reg1_int_udsmean,
                 reg2_int_udsmean,
                 reg1_int2_udsmean,
                 reg2_int2_udsmean)
se_uds <- list(coef_reg1_int_udsmean[,2],
               coef_reg2_int_udsmean[,2],
               coef_reg1_int2_udsmean[,2],
               coef_reg2_int2_udsmean[,2])
p_uds <- list(coef_reg1_int_udsmean[,4],
              coef_reg2_int_udsmean[,4],
              coef_reg1_int2_udsmean[,4],
              coef_reg2_int2_udsmean[,4])
covlabs_uds <- c("UDS democracy ($t-5$)",
                 "OECD member ($t-5$)",
                 "GDP/capita, \\$ (ln, $t-5$)",
                 "Population (ln, $t-5$)",
                 "Trade, \\% GDP ($t-5$)",
                 "Urban population, \\% ($t-5$)",
                 "Kyoto ratification ($t-5$)",
                 "EU member ($t-5$)",
                 "UDS $\\times$ OECD member ($t-5$)",
                 "UDS $\\times$ GDP/capita (ln) ($t-5$)")
n_uds <- list(c("Countries","70","68","68","68"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_uds,
          header = F,
          list_uds,
          covariate.labels = covlabs_uds,
          omit.stat = omit.stats,
          add.lines = n_uds,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > India/China-excluded regression ####

reg1_int_minusindiachina <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5,
                                   data = coal[coal$country != "India" & coal$country != "China",], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_minusindiachina <- coeftest(reg1_int_minusindiachina, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_minusindiachina <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5 +
                                     log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                                     urban_pop_l5 +
                                     kyoto_l5 + EU_l5,
                                   data = coal[coal$country != "India" & coal$country != "China",], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_minusindiachina <- coeftest(reg2_int_minusindiachina, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_minusindiachina <- plm(log(cumul_capacity_post80) ~ polity2_l5*log(gdp_pc_l5),
                                    data = coal[coal$country != "India" & coal$country != "China",], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_minusindiachina <- coeftest(reg1_int2_minusindiachina, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_minusindiachina <- plm(log(cumul_capacity_post80) ~ polity2_l5*log(gdp_pc_l5) +
                                      log(population_l5) + trade_pctgdp_l5 +
                                      urban_pop_l5 +
                                      kyoto_l5 + EU_l5 + oeec_oecd_l5,
                                    data = coal[coal$country != "India" & coal$country != "China",], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_minusindiachina <- coeftest(reg2_int2_minusindiachina, vcov=function(x) vcovHC(x, cluster = "group")) )

list_indchn <- list(reg1_int_minusindiachina,
                    reg2_int_minusindiachina,
                    reg1_int2_minusindiachina,
                    reg2_int2_minusindiachina)
se_indchn <- list(coef_reg1_int_minusindiachina[,2],
                  coef_reg2_int_minusindiachina[,2],
                  coef_reg1_int2_minusindiachina[,2],
                  coef_reg2_int2_minusindiachina[,2])
p_indchn <- list(coef_reg1_int_minusindiachina[,4],
                 coef_reg2_int_minusindiachina[,4],
                 coef_reg1_int_minusindiachina[,4],
                 coef_reg2_int2_minusindiachina[,4])
covlabs_indchn <- c("Polity2 ($t-5$)",
                    "OECD member ($t-5$)",
                    "GDP/capita, \\$ (ln, $t-5$)",
                    "Population (ln, $t-5$)",
                    "Trade, \\% GDP ($t-5$)",
                    "Urban population, \\% ($t-5$)",
                    "Kyoto ratification ($t-5$)",
                    "EU member ($t-5$)",
                    "Polity2 $\\times$ OECD member ($t-5$)",
                    "Polity2 $\\times$ GDP/capita (ln) ($t-5$)")
n_indchn <- list(c("Countries","69","67","67","67"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_indchn,
          header = F,
          list_indchn,
          covariate.labels = covlabs_indchn,
          omit.stat = omit.stats,
          add.lines = n_indchn,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > 1990s-excluded regression ####

reg1_int_minus90s <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5,
                            data = coal[coal$year < 1990 | coal$year > 1999,], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_minus90s <- coeftest(reg1_int_minus90s, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_minus90s <- plm(log(cumul_capacity_post80) ~ polity2_l5*oeec_oecd_l5 +
                              log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                              urban_pop_l5 +
                              kyoto_l5 + EU_l5,
                            data = coal[coal$year < 1990 | coal$year > 1999,], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_minus90s <- coeftest(reg2_int_minus90s, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int2_minus90s <- plm(log(cumul_capacity_post80) ~ polity2_l5*log(gdp_pc_l5),
                             data = coal[coal$year < 1990 | coal$year > 1999,], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int2_minus90s <- coeftest(reg1_int2_minus90s, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int2_minus90s <- plm(log(cumul_capacity_post80) ~ polity2_l5*log(gdp_pc_l5) +
                               log(population_l5) + trade_pctgdp_l5 +
                               urban_pop_l5 +
                               kyoto_l5 + EU_l5 + oeec_oecd_l5,
                             data = coal[coal$year < 1990 | coal$year > 1999,], model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int2_minus90s <- coeftest(reg2_int2_minus90s, vcov=function(x) vcovHC(x, cluster = "group")) )

list_minus90s <- list(reg1_int_minus90s,
                      reg2_int_minus90s,
                      reg1_int2_minus90s,
                      reg2_int2_minus90s)
se_minus90s <- list(coef_reg1_int_minus90s[,2],
                    coef_reg2_int_minus90s[,2],
                    coef_reg1_int2_minus90s[,2],
                    coef_reg2_int2_minus90s[,2])
p_minus90s <- list(coef_reg1_int_minus90s[,4],
                   coef_reg2_int_minus90s[,4],
                   coef_reg1_int2_minus90s[,4],
                   coef_reg2_int2_minus90s[,4])
covlabs_minus90s <- c("Polity2 ($t-5$)",
                      "OECD member ($t-5$)",
                      "GDP/capita, \\$ (ln, $t-5$)",
                      "Population (ln, $t-5$)",
                      "Trade, \\% GDP ($t-5$)",
                      "Urban population, \\% ($t-5$)",
                      "Kyoto ratification ($t-5$)",
                      "EU member ($t-5$)",
                      "Polity2 $\\times$ OECD member ($t-5$)",
                      "Polity2 $\\times$ GDP/capita (ln) ($t-5$)")
n_minus90s <- list(c("Countries","70","68","68","68"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_minus90s,
          header = F,
          list_minus90s,
          covariate.labels = covlabs_minus90s,
          omit.stat = omit.stats,
          add.lines = n_minus90s,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity in Megawatts (ln)")
)

# > Unified Democracy Scores simulations: 1000 draws from UDS posterior, interaction with GDPpc ####

rep.col <- function(x, n){
  matrix(rep(x,each=n), ncol=n, byrow=TRUE)
}
.uds_draw_results <- rep.col(1:18,1000)
.uds_draw_results <- as.data.frame(.uds_draw_results)
.uds_draw_results[.uds_draw_results > 0] <- NA

for (i in 1:1000){
  .uds_draw <- coeftest(lm(log(cumul_capacity_post80) ~ as.matrix(coal[21+i])*log(gdp_pc_l5) +
                             log(population_l5) + trade_pctgdp_l5 +
                             urban_pop_l5 +
                             kyoto_l5 + EU_l5 + oeec_oecd_l5 + factor(country) + factor(year),
                           data = coal), vcov=function(x) vcovHC(x, cluster = "group"))
  .uds_draw_results[1:9,i] <- .uds_draw[c(2:9,102),1]
  .uds_draw_results[10:18,i] <- .uds_draw[c(2:9,102),4]
}

uds_draw_results <- as.data.frame(t(.uds_draw_results))

uds_rs_results <- data.frame(mean=rep(NA,9), a05=rep(NA,9), a01=rep(NA,9))

for (i in 10:18){
  uds_rs_results$mean[i-9] <- mean(uds_draw_results[,i-9])
  uds_rs_results$a05[i-9] <- mean(uds_draw_results[,i] <= .05)
  uds_rs_results$a01[i-9] <- mean(uds_draw_results[,i] <= .01)
}

rownames(uds_rs_results) <- c("UDS (t-5)","GDP/capita, $ (ln, t-5)","Population (ln, t-5)",
                              "Trade, % GDP (t-5)", "Urban population, % (t-5)","Kyoto ratification (t-5)",
                              "EU member (t-5)","OECD member (t-5)","UDS * GDP/capita (t-5)")
colnames(uds_rs_results) <- c("Mean coefficient", "% p < .05", "% p < .01")

print(xtable(uds_rs_results),floating=FALSE,latex.environments="center")

# > Unified Democracy Scores simulations: 1000 draws from UDS posterior, interaction with OECD ####

.uds_draw_results2 <- rep.col(1:18,1000)
.uds_draw_results2 <- as.data.frame(.uds_draw_results2)
.uds_draw_results2[.uds_draw_results2 > 0] <- NA

for (i in 1:1000){
  .uds_draw <- coeftest(lm(log(cumul_capacity_post80) ~ as.matrix(coal[21+i])*oeec_oecd_l5 + log(gdp_pc_l5) +
                             log(population_l5) + trade_pctgdp_l5 +
                             urban_pop_l5 +
                             kyoto_l5 + EU_l5 + factor(country) + factor(year),
                           data = coal), vcov=function(x) vcovHC(x, cluster = "group"))
  .uds_draw_results2[1:9,i] <- .uds_draw[c(2:9,102),1]
  .uds_draw_results2[10:18,i] <- .uds_draw[c(2:9,102),4]
}

uds_draw_results2 <- as.data.frame(t(.uds_draw_results2))

uds_rs_results2 <- data.frame(mean=rep(NA,9), a05=rep(NA,9), a01=rep(NA,9))

for (i in 10:18){
  uds_rs_results2$mean[i-9] <- mean(uds_draw_results2[,i-9])
  uds_rs_results2$a05[i-9] <- mean(uds_draw_results2[,i] <= .05)
  uds_rs_results2$a01[i-9] <- mean(uds_draw_results2[,i] <= .01)
}

rownames(uds_rs_results2) <- c("UDS (t-5)","OECD member (t-5)","GDP/capita, $ (ln, t-5)","Population (ln, t-5)",
                               "Trade, % GDP (t-5)", "Urban population, % (t-5)","Kyoto ratification (t-5)",
                               "EU member (t-5)","UDS * OECD (t-5)")
colnames(uds_rs_results2) <- c("Mean coefficient", "% p < .05", "% p < .01")

print(xtable(uds_rs_results2),floating=FALSE,latex.environments="center")

# > Subtracting decommissioned capacity ####

reg1_int_no.decom <- plm(log1p(cumul_capacity_post80_nodecom) ~ polity2_l5*oeec_oecd_l5,
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_no.decom <- coeftest(reg1_int_no.decom, vcov=function(x) vcovHC(x, cluster = "group")) )

reg1_int_no.decom2 <- plm(log1p(cumul_capacity_post80_nodecom) ~ polity2_l5*log(gdp_pc_l5),
                             data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg1_int_no.decom2 <- coeftest(reg1_int_no.decom2, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_no.decom <- plm(log1p(cumul_capacity_post80_nodecom) ~ polity2_l5*oeec_oecd_l5 +
                              log(gdp_pc_l5) + log(population_l5) + trade_pctgdp_l5 +
                              urban_pop_l5 +
                              kyoto_l5 + EU_l5,
                            data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_no.decom <- coeftest(reg2_int_no.decom, vcov=function(x) vcovHC(x, cluster = "group")) )

reg2_int_no.decom2 <- plm(log1p(cumul_capacity_post80_nodecom) ~ polity2_l5*log(gdp_pc_l5) +
                               log(population_l5) + trade_pctgdp_l5 +
                               urban_pop_l5 +
                               kyoto_l5 + EU_l5 + oeec_oecd_l5,
                             data = coal, model = "within", effect = "twoways", index = c("country", "year"))
( coef_reg2_int_no.decom2 <- coeftest(reg2_int_no.decom2, vcov=function(x) vcovHC(x, cluster = "group")) )

list_nodecom <- list(reg1_int_no.decom,
                     reg1_int_no.decom2,
                     reg2_int_no.decom,
                     reg2_int_no.decom2)
se_nodecom <- list(coef_reg1_int_no.decom[,2],
                   coef_reg1_int_no.decom2[,2],
                   coef_reg2_int_no.decom[,2],
                   coef_reg2_int_no.decom2[,2])
p_nodecom <- list(coef_reg1_int_no.decom[,4],
                  coef_reg1_int_no.decom2[,4],
                  coef_reg2_int_no.decom[,4],
                  coef_reg2_int_no.decom2[,4])
covlabs_nodecom <- c("Polity2 ($t-5$)",
                     "OECD member ($t-5$)",
                     "Polity2 $\\times$ OECD member ($t-5$)",
                     "GDP/capita, \\$ (ln, $t-5$)",
                     "Polity2 $\\times$ GDP/capita (ln) ($t-5$)",
                     "Population (ln, $t-5$)",
                     "Trade, \\% GDP ($t-5$)",
                     "Urban population, \\% ($t-5$)",
                     "Kyoto ratification ($t-5$)",
                     "EU member ($t-5$)")
n_nodecom <- list(c("Countries","71","69","69","69"))

stargazer(title = "Democracy and Commissioned Coal Plant Capacity",
          se = se_nodecom,
          header = F,
          list_nodecom,
          covariate.labels = covlabs_nodecom,
          omit.stat = omit.stats,
          add.lines = n_nodecom,
          no.space = T,
          suppress.errors = F,float = F,
          dep.var.labels = c("Cumulative Commissioned Capacity, Less Decommissioned Capacity, in Megawatts (ln)")
)
